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ABSTRACT 

The velocity dispersion, or eccentricity distribution, of protoplanets interacting with planetesimals 
is set by a balance between dynamical friction and viscous stirring. We calculate analytically the 
eccentricity distribution function of protoplanets embedded in a cold, shear-dominated planetesimal 
swarm. We find a distinctly non-Rayleigh distribution with a simple analytical form. The peak of 
the distribution lies much lower than the root-mean-squared value, indicating that while most of the 
bodies have similarly small eccentricities, a small subset of the population contains most of the thermal 
energy. We also measure the shear-dominated eccentricity distribution using numerical simulations. 
The numerical code treats each protoplanet explicitly and adds an additional force term to each 
body to represent the dynamical friction of the planetesimals. Without fitting any parameters, the 
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eccentricity distribution of protoplanets in the N-body simulation agrees with the analytical results. 
This distribution function provides a useful tool for testing hybrid numerical simulations of late-stage 
planet formation. 
Subject headings: planets and satellites: formation — solar system: formation 



1. INTRODUCTION 



> 

00 

Terrestrial planets, ice giants, and the cores of the gas giants are thought to form by accretion of planetesimals into 
r^f- protoplanets. The protoplanets emerge from the swarm of planetesimals after an epoch of runaway accretion. The 
f^ | subsequent dynamics of the protoplanets set several important features of the final planetary configuration, such as 
VO . the mass and number of planets or cores. It is difficult to constrain this evolution, however, without constraining the 
^D ' properties of the disk in which they are embedded. 

One important yet uncertain parameter is the size of the planetesimals, the building blocks. The outer Solar System 
and the later stages of formation in the inner Solar System likely lack gas, allowing the formation of kilometer-size 
bodies through gravitational instabilities. Those bodies collide and grind each other down to even smaller sizes in 
a collisional cascade. The existence of bodies small enough to damp their own velocity dispersion is an inevitable 
£3 . conclusion from the existence of Uranus and Neptune (Goldrcich ct al. 2004b). Without such small bodies, the ability 
of a growing protoplanet to gravitationally focus the planetesimals becomes inefficient, and the growth timescale 



becomes too long, of order 10 12 years in the outer solar system. 

The unavoidable influence of the planetesimals make numerical studies of planet formation difficult to carry out 
accurately. Despite modern computational power, an integration of the equations of motion for each body in a 
protoplanet and planetesimal swarm is impossible. Even without allowing planetesimal frag mentation, the n umber 
of kilometer-size bodies needed to comprise a Neptune size mass is humongous, of order 10 12 . Kokubo & Ida (1996) 
performed numerically feasible but physically less appropriate N-body simulations of a proto-planetary disk in which 
the size of planetesimals is larger than the value required to form the ice giants of our Solar System. Although 
interesting from a dynamical viewpoint, the results of such simulations can not be extrapolated to the scenario of 
smaller planetesimals since they lack collisional damping. 

An alternative numerical approach to studying these systems is a coagulation code l|Leeil2000i iKenvon fc Luul ll998) 
in which the bodies are divided into size bins and the interaction of each pair of bins is calculated statistically. This 
approach fails once the number of bodies in any bin is not sufficiently large. IKenvon fc Brom lcv (2006) have developed 
a hybrid code that treats planetesimals statistically while a small number of large bodies are integrated individually. 

In this paper, we examine the processes that shape the eccentricity distribution of the large bodies. We assume, 
simply, that the planetesimals constitute a cold disk due to sufficiently frequent collisions. As a first step, we include 
the dynamical friction that the planetesimals exert on the large bodies but ignore the much slower process of their 
accretion onto those bodies. The rates of cooling from dynamical friction and heating from mutual excitations are 
discussed in |2] We write a Boltzmann equation to show the change in the distribution function of eccentricities due 
to each process in SJ3 and discuss the solution to that equation. In 0we present the results of complementary N-body 
simulations designed to measure the eccentricity distribution directly. A discussion of the results follows in 50 

2. SHEAR-DOMINATED COOLING AND HEATING RATES 

The eccentricities of the protoplanets represent a kind of "thermal" energy in their orbits, relative to perfectly 
circular motion. The extra non-circular velocity itself varies in magnitude and direction over an orbital period; it is 
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simpler to use the eccentricity, a constant of motion for the two-body problem. Specifically, we calculate the vector 
eccentricity, 

vxH r . , 

GM P r 

This expression relates the eccentricity of the particle, e, to the particle's position, r, its velocity, v, its orbital angular 
momentum vector, H, and its mass, M p . In general, a protoplanet can have an inclination relative to the disk plane, 
and the eccentricity vector can have three components. However, we show in H2. 41 that the shear-dominated regime 
strongly inhibits the growth of inclinations. Two dimensions then suffice to describe the configuration space of e. 
We use the quantity of the Hill radius repeatedly in this work; for reference we define its value as 

/ M \ 1/3 
Rh = {^) a = R/a, (2) 

where M p is the mass of a particle, a is its semi-major axis, R is its radius, and 

/ Mp \ R I Pp \ R 



The Hill radius in turn specifies an eccentricity, the Hill eccentricity, 

e H = Rh/o,. (4) 

We restrict this study to disks where the majority of the bodies have eccentricities lower than e#, known as the 
shear-dominated regime. ^^^^^^_^^^^^^_ ^_^ 

For most of this paper, we employ the "two groups" approximation IjWetherill fc StewartHl989t iGoldreich et alJ 
2004b) and split the disk into two uniform populations. One group is the numerous smaller bodies, or "planetesimals." 
We denote their surface mass density as a. The other group, the "protoplanets" , consists of the bodies that dominate 
the excitations of the disk particles. Each protoplanet has a radius -R, mass M, and eccentricity e. We write the total 
surface mass dens ity in protoplanets as S. We assume a > S although this is not always true in the final stages of 
planet formation l|Goldreich et al.H2004a|) . 

2.1. Eccentricity Excitation of Protoplanets 

We analyze the interaction of two protoplanets from a frame rotating with a reference orbit at a semi-major axis a. 
The difference between the Keplerian angular velocity at each radius induces a shearing motion between particles on 
nearby circular orbits. For an orbit interior to a by a distance b, 

Q 7 

fl Tel (b) = Cl(a + b) -0(a) « -0-, (5) 

2 a 

in the limit of b <C a. This angular frequency also specifies the rate of conjunctions for the two bodies with orbits 
separated by b. 

The change in their eccentricity from each conjunction can be calculated analytically for two nearly circular orbits 
when b 3> Rh'- 

e k = A k e H ( -5— ) , (6) 
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Kg and Jfi are mo dified Bessel functions of the second kind l|Goldreic h fc Tremainei 119801 iPetit & Henonl Il98fit 
iDuncan et alJll989l) . We note that e k refers to the perturbed body, while en and Rh refer to the Hill parame- 
ters of the perturbing one. The kick, viewed as a change in the eccentricity vector, is independent of the original 
eccentricity of the particle. Its orientation is perpendicular to the line connecting the two protoplanets and the sun at 
conjunction; we therefore assume it is random. 

The eccentricity kick given by one protoplanet is strongest for the particles that approach with an impact parameter 
on the order of the Hill radius. Interactions from a greater distance, however, occur more often. In a shear-dominated 
disk, eccentricities are small, e <§; en ; to change that eccentricity significantly only requires sm all perturbations. These 
frequent but weaker perturbations dictate the overall velocity evolution of the protoplanets (Goldrei ch et alJ|2004bt 
Rafikov 2004). 

Specifically, the average differential rate that one protoplanet receives eccentricity kicks of strength e from other 
protoplanets is given by 

d^ox(e) = 2 n big -06(e) — de, (8) 



where ribi g is the number surface density of protoplanets and |fi6(e) is the velocity of encounters at those separations 
(given by eq. |3J). The factor of two accounts for the combination of interior and exterior encounters. The excitation 
rate of a protoplanet with eccentricity e is then the rate of kicks comparable in magnitude to its current eccentricity: 
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The inverse of this rate can be interpreted as the timescale for a protoplanet 's eccentricity to change by an amount e. 

2.2. Dynamical Friction 

As each protoplanet moves through the disk, it scatters and excites the eccentricities of the planetesimals that 
surround it. Cold planetesimals that approach a protoplanet with impact parameters of about a Hill sphere leave 
with ~ m en of additional momentum. This can either add to or subtract from the eccentricity of the protoplanet 
depending on the relati ve orientation between the pre-encounter eccentricities of the protoplanet and planetesimal. 
We write the net effect l)Goldreich et al.ll2004bT) 
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-n s R H eHm(e H + e) + n s R H e H m{e H - e), (10) 

for a number surface density of planetesimals n s . This formula yields the damping rate, or the inverse damping time, 
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Calculating the coefficie nt Cd requires a more precise analysis of planetesimal scattering. We adopt the value Cd 
found bv lOhtsuki et afl l)2002() . who measure the coefficient numerically. 

2.3. Planetesimal Interactions 

The distribution of planetesimal eccentricities does not affect our results, as neither the excitation nor the damping 
rates depend on their eccentricity as long as the planetesimals remain in a shear-dominated state. In this work, we 
focus on a range of parameters such that collisional cooling keeps the eccentricities of planetesimals below eu and 
enforces the condition of shear domination. 

2.4. Inclinations 

An orbit with a small inclination angle i carries its particle out of the disk plane on vertical excursions of size ~ ia. 
An interaction that excites that particle's eccentricity also affects its inclination, but with a magnitude inhibited by 
the geometry of the distant encounter: 
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where b is the impact parameter of the perturber and i^ is the resulting change in inclination from an encounter. In 
contrast, planetesimals just entering the Hill sphere of a protoplanet damp the protoplanet's non-circular velocity in all 
dimensions; no equivalent geometric factor inhibits the damping of inclinations. Wit h the growt h of inclinatio ns sup- 
pressed, shear-domin ated protoplanet disks are effectively two-dimensional ( Wctherill & Stewart 1993; Goldreich et al. 
l2004btlRalil^l20T)l . 

2.5. The Eccentricity Distribution - a Qualitative Discussion 

The dynamical friction rate sets a characteristic time over which the eccentricities of all of the bodies are changed 
significantly. In this sense, the eccentricity distribution of the proto-planetary swarm is reset every t<j. The excitation 
rate, however, varies with e. Equating the excitation rate, equation |5J and the damping rate, equation II II yields an 
important reference value, e eq : 
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Statistically, each protoplanet receives one kick of this magnitude every damping timescale. 

We deduce the distribution of eccentricities on each side of e oq by examining the dependence of the kicking rate on 
eccentricity. Excitation to e ^ e cq requires a kick e/. 3> e cq . Such strong kicks occur less often in one damping timescale 
than kicks of strength e oq by a factor of e cq /efc. With fewer kicks to populate the high eccentricity distribution, the 
numbe r of bodies with such eccentricities echoes the rate of kicks and falls off with eccentricity as e _1 ( Goldreich ct al. 
20043). 

Kicks of order e <C e eq occur frequently in each damping timescale, thereby overwhelming the effects of dynamical 
friction on the lowest eccentricity bodies. A dynamic equilibrium dominated by only the stirring mechanism implies 
that kicks to and from every eccentricity vector occur at the same rate. For this to be true the distribution must be 
constant over the configuration space. The number of bodies with an eccentricity of order e <C e cq then scales as the 
area of configuration space available, ~ e 2 . 



3. A BOLTZMANN EQUATION 

In the following section we develop a differential equation to describe analytically the distribution function of 
protoplanet eccentricities. We construct this equation in the spirit of the Boltzmann equation, examining the change 
in the number of bodies with a particular eccentricity due to the effects of dynamical friction and viscous stirring. 

The space of possible eccentricities is inherently two-dimensional (eq. [IJ, since inclinations can be neglected ({j 
Additionally the interaction rates depend only on the magnitude of the protoplanet eccentricity forcing the distribution 
function to share this dependence: /(e) = /(e). The two-dimensional /(e) is related to the number of bodies with 
velocity on the order of e by its integral, roughly e 2 /(e). 

Dynamical friction lowers the eccentricities of all bodies proportionally to their eccentricity. Equivalently, the number 
of bodies with a certain e changes as the protoplanets with that value are damped to lower eccentricities and replaced 
by bodies from a higher eccentricity. We write this as: 
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(14) 



where we have used de/dt = —e/rd for the effects of dynamical friction. 

At a given e, particles are kicked to a new eccentricity e„ at an average rate that depends on the magnitude of the 
kick, |e„ — e|. Also, particles with an initial eccentricity e„ are kicked to e at the same rate. The total flux of particles 
to and from a given eccentricity is: 
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where p(e) describes the rate at which bodies experience changes in their eccentricities by an amount e. This is the 
two-dimensional analog of the excitation rate equation [5J 
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The sum of the dynamical friction terms and the kicking integral describes the dynamics of shear-dominated proto- 
planets interacting with each other in a smooth disk of planetesimals. The combined influence of these two processes 
can bring the protoplanets into an equilibrium state, where the number of particles with eccentricity e remains constant 
in time: 



= ^P- + ^ + / / P(|e„ - e|) [/(e„) - /(e)] d 2 e r , 
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3.1. The Solution 



We show in Appendix A that 
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This function is the equilibrium eccentricity distribution of 
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satisfies the equilibrium equation, equation 1171 for all e. 
shear-dominated protoplanets. 

The solid line in Fig. 1 shows a distribution function for E w 0.002 g cm~ 2 and a = 0.1 g cm~ 2 . Although the 
function formally extends above e#, we stress that it is only accurate for eccentricities e <C en- Both the dynamical 
friction and the excitation rates feas. 1141 &: I16|) are not valid for e > e#. 

Several moments of the distribution can be calculated in terms of the only free parameters: 
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According to equation 1181 (e) is infinite. However, the largest single kick in eccentricity from an almost circular 
protoplanet encounter is of order en ■ Truncating the integral at e# produces the logarithmic term in the expression 
above. Moments higher than the mean also diverge; realistically, they are dominated by the bodies with the highest 
eccentricities, of order en- 

It is easy to see that this solution, in the high- and low-eccentricity limits, produces the same power-laws discussed 
in H2.5I In fact, it can be shown directly from equation 1171 that any solution to this differential equation reduces to 
those limiting power-laws. 



4. NUMERICAL SIMULATIONS 

Here we describe a direct measurement of the eccentricity distribution from gravitational N-body simulations that 
include an additional force to represent dynamical friction. 

The N-body part of our simulation uses Gauss's equations to evolve a set of orbital constants chosen to vary slowly 
under small perturbations. A modified version of Kepler's equation produces the orbital ph ase for each body at each 
time step. The IDA solver from the SUNDIALS software package (Hindma rsh et al]l2005[) integrates the dynamical 
equations. During close encounters of two protoplanets, we integrate their motion relative to the center of mass of the 
pair. 

We represent the planetesimal population in these simulations with an extra force term that damps the non-circular 
velocities of the protoplanets at the rate lZd{e), given by eauationlTTl An ad-hoc transition between the damping rate 
for e < en and the appropriate rate for e > en prevents unphysical enhancements of the damping force during close 
encounters. The growth of protoplanets in mass due to planetesimal accretion is not incl uded; the accretion rate is 
always lower than the dynamical friction rate and will not affect the eccentricity evolution CGol dreich et alJl200 4b1 . 

Each simulation begins with the protoplanets on circular orbits with random phases and random semi-major axes, 
within a chosen annulus. The average spacing between bodies, M/(£27ra), is several Hill radii. The protoplanets 
interact for several damping timescales Td before the distribution reaches equilibrium. 

We record the eccentricity of the protoplanets every At w OTr^ starting at about 100 Td- Several hundred orbits 
produces a well-populated histogram of eccentricities. The bodies in the inner and outer edges of the disk are not 
measured, to avoid artificial boundary effects that inhibit excitations. We bin the resulting eccentricities logarithmi- 
cally. Errors are assigned to each bin according to a Poisson distribution with the sample size defined as the product 
of the number of bodies measured and the sampling time in units of the damping timescale Td- Since each protoplanct 
suffers a significant change in eccentricity every Td, one measurement of the eccentricity distribution is independent 
from a previous measurement if they are separated in time by Td- We sample faster than Td to increase the resolution 
of the histogram slightly. 

The statistical error bars do not take into account the inhomogeneity of the protoplanet disk on small length scales. 
Given a surface density, the mass of a single protoplanet sets a typical radial separation between bodies. This length 
scale corresponds to an eccentricity scale through equation [S] (in the simulations presented here, this value is slightly 
below en ). As the disk evolves, the viscous stirring causes migrations in the semi-major axes of the particles that 
smooth the average radial distribution. If measured only over intervals shorter than the migration timescale, the 
eccentricity distribution may vary for eccentricities above the eccentricity set by the typical separation. Fluctuations 
from this effect are visible in Figs. 1 and 2. 

Several simulations of disks with different protoplanct mass distributions are presented below. 

4.1. Equal Mass Protoplanets 

Figure 1 shows the eccentricity distribution measured from a simulated disk of 120 equal mass protoplanets (M = 
2.5 x 10 -9 Mq) with surface densities £ sw 0.002 g cm~ 2 and a = 0.1 g cm -2 . A single population of protoplanets best 
reflects the "two groups" approximation we use to derive equations El and 1111 The analytic solution, equation 1181 for 
the same parameters in the simulation is superposed on Figure 1. While the overall match is not perfect, the shape 
of each curve is strikingly similar. The two curves match extremely well if one is shifted by around 15 percent in the 
e direction. This difference is attributable to the difficulty of assigning a correct value of S to the simulation given a 
finite number of protoplanets. 

We note that there are no free parameters in this comparison. The numerical distribution is a direct counting of the 
number of bodies within each eccentricity bin, while a choice of S, <j, and M completely specifies the analytical curve. 

4.2. Mass Distributions 

Naturally occurring protoplanet populations exhibit non-trivial distributions in mass. Before describing such a disk 
in the framework we have developed, we clarify several points. 

Protoplanets with different masses, or equivalently, different radii, experience different viscous stirring rates. We 
decompose the total surface density in protoplanets, S, into a differential quantity, dT,/dR, and write the excitation 
rate of a body with radius R as 

[(Efll e H (R>) 

n * {e > R) ~JdR7p-E7^^ dR - (20) 

The identity ejj(R') — (R 1 /R)eu{R) when substituted into equation 1201 yields 

Q 1 e H (R) f dS 

^^-W^eJW^- (21) 

In words, the excitation rate of one body only depends on t he total surface d ensity of all other bodies, regardless of the 
specific mass distribution. This differs from the assertion by Goldrcich et al. (2004b) that only the most massive bodies 
contribute to the viscous stirring rate. Equation 12 II seems to indicate that there should be no distinction between big 
bodies and small bodies since every body contributes to the viscous stirring. A closer investigation uncovers the mass 
range of bodies that provide significant stirring. 
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FIGURE 1: A plot of equation 1181 superposed with the results of a numerical simulation. The simulated disk contains 120 bodies of 
mass M = 5 X 10 24 g, orS fs 0.002 g cm -2 . A planetesimal surface density of a = 0.1 g cm -2 is included. We assume each bin obeys 
Poisson statistics and assign errors based on a population size of N b N Td , where N b is the number of bodies in the simula tion , and N Td 
is the duration of the simulation in units of damping timescales. The solid line shows the distribution as given by equation 1181 using the 
same values of £ and <r. A Rayleigh distribution with a similar peak eccentricity is plotted with the dashed line. 



Eccentricity kicks of strength e^ can occur at any combination of M and b that satisfies the inverse square law of 
gravitation: e^ ~ M(R')b(R')~ 2 . However, the smallest impact parameter that contributes to a body's excitation is 
about Rh . A minimum b(R') ~ Ru sets a minimum mass for bodies to kick a body with mass M by an amount e^: 



A/ mi n(e fe ,i?)~— M 
Likewise, a body can only be as far away as its radial position in the disk, a. This sets a maximum mass, 
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For a choice of the most relevant kick strength, e^, these limits define the sizes of bodies that participate in the 
excitation of a body with size R. 

As a numerical confirmation of these results, we simulate a disk of planetesimals with a surface mass density 
g = 0.2 g cm~ 2 and 120 protoplanets. In this case, we divide the protoplanets into two groups of different mass: sixty 
of mass mi — 2x 10 24 g, and sixty of mass m,2 — 3.8 x 10 25 g. These masses are within the limits set by equations I22I&: 
1231 We plot the absolute eccentricity distribution of each mass group binned separately in Figure 2. Additionally, we 
plot the analytic distributions given by a, as specified above, and S, the sum of the surface densities of both groups. 

It is clear that each group of protoplanets with the same mass matches the analytic distribution well. The offset 
between the peak of each group is due to the dependence of the distribution on the Hill eccentricity of each body. 
In general, the distribution for a swarm of protoplanets with a mass distribution is merely the sum of individual 
distributions for protoplanets of each mass. 

5. CONCLUSIONS 

We presented an analytic model for the distribution function of the eccentricities of a protoplanet population em- 
bedded in a shear-dominated planetesimal disk. The eccentricity distribution measured with numerical simulations 
matches the analytic result very well. 

Since we have manually inserted the dynamical friction rate that we expect into the numerical simulations, this work 
does not test our prescription of dynamical friction. However, the numerical and analytic representations of viscous 
stirring are completely independent. Eauation ll7l uses a viscous stirring rate involving only distant encounters. In our 
numerical simulations, Newton's laws dictate the protoplanet interactions directly without any simplifying assumptions. 
The consistency of the two calculations proves that in a two-dimensional shear-dominated disk, interactions between 
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FIGURE 2: A comparison of the results of a numerical simulation of protoplanets in a perfectly bimodal mass distribution (mi = 
2 X 10 24 g, m,2 = 3.8 X 10 25 g). We simulate sixty bodies of each mass, for a total surface density in protoplanets of S m 0.003 g cm -2 and 
a planctcsimal sur face density of a = 0.2 g cm -2 . The eccentricities of each mass group are binned separately; each distribution is a good 
match to equation 1181 when scaled to the appropriate Hill eccentricity. The error bars are assigned following the same algorithm as Figure 
1. 



non-crossing orbits are entirely responsible for setting the eccentricities of the protoplanets. The analytic form of the 
distribution function provides a simple way to test similar hybrid numerical codes. 

Several features of the distribution highlight interesting properties of the dynamics. We reason in H2. 51 that most 
protoplanets have eccentricities ~ (T,/a)eH, the value of e where the excitation and damping timescales are equal. The 
distribution function shows this to be true: the median and mean (up to a logarithmic factor) of any distribution are on 
the order of this equilibrium eccentricity. Higher moments of the distribution, however, are dominated by the highest 
eccentricity bodies. This signals that different statistics of the distribution can reflect different subsets of the overall 
population. For example, the average "thermal" energy of the protoplanets is represented by the root-mean-squared 
eccentricity, (e 2 ). The fractionally fewer bodies with eccentricities close to e# dominate (e 2 ) and thus, contain most 
of the energy. 

The shape of the distribution also merits discus sion. N-body integra tions of a group of single mass bodies show that 
their eccentricities follow a Rayleigh distribution ((Ida fc Makinolll992lh For reference, we plot a Rayleigh distribution 
in Figure 1. It is entirely inconsistent with our calculations. This is not surprisin g. In addit i on to simulating bodies in 
the regime of eccentricities that are large compared to the Hill eccentricity, Ida & Makino (1992) do not include any 
effects that can balance the mutual excitations of their particles. The dynamical friction in our simulations balances 
the viscous stirring and establishes the equilibrium distribution we derive. 



RS is an Alfred P. Sloan Fellow and a Packard Fellow. 
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APPENDIX 

THE ANALYTIC DISTRIBUTION FUNCTION 

Here we outline the evaluation of the equilibrium equation, equation 1171 using the distribution function, equation 
1181 To simplify the notation, we rescale all eccentricities by e* and algebraically manipulate the coefficients of each 
term in eauation ll7l We are left with the equivalent burden of proving that 



3 (e) = (l + e 2 )- 3/2 
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The left-hand side is easy to compute: 
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(Al) 
(A2) 

(A3) 



To integrate of the right-hand side, we translate the origin of the integration variables by e and rotate them to align 
e with one of the coordinate axes. In those coordinates: 



/ 



oo />oo 



i 



oc J — oo 



(k 2 + ^2)3/2 



1 



(1 + e 2 ) 3 /2 (1 + k 2 + {h + e) 2 ) 3 / 2 



dk dh, 



(A4) 



with e n = {k. h}. 
After the integration over k, we rewrite the integral in terms of a new variable h! = (1 4- e 2 )/(e/i), 



I = 



2e 



(l + e 2)5/2 (l + e 2) 



s ^(eVd + e!)) WUh , 



dz 2 



(A5) 



where, z = —2h! — h' 2 , and E{e 2 z/(1 + e 2 )) is the complete elliptic integral of the second kind. 

We change the integration variable to z, taking care to evaluate the integrand with the appropriate branch of the 
double- valued relation h'(z). The integral evaluates to 



/ = 



-iii 



(l + e 2 )3 (i + e 2)2y o L^/T 



4-3z 



d 2 E(e 2 z/(l + e 2 )) 
dz 2 



dz . 



(A6) 



With the second derivative of the elliptic function expressed as a power series, each term can be integrated over z. 
The remaining power series in e 2 /(l + e 2 ) equals 



I = 



4tt 



37re 4 



(1 



o2^3 



(1 



,2U 



2-ri 



Fi -,1;3 



) -"-} "I 



(1 



"2 2Fl (2' 2;3; (TT^) 



(A7) 



After additional algebraic manipulation, this result equals the left hand side of the original equation (eq. IA3 



